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We have developed an improved version of the quantum transfer matrix algorithm. The extreme 
eigenvalues and eigenvectors of the transfer matrix are calculated by the recently developed look- 
ahead Lanczos algorithm for non-Hermitian matrices with higher efficiency and accuracy than by 
the power method. We have applied this method to the Heisenberg ladder. The temperature 
dependence of the susceptibility, specific heat, correlation length and nuclear spin relaxation rate 
1/Ti are calculated. Our results support the existence of a spin gap of about 0.5 J. 



I. INTRODUCTION 

The behavior of one-dimensional (ID) strongly corre- 
lated systems and spin chains is by now quite well un- 
derstood. For two-dimensional (2D) strongly correlated 
systems there are many open questions. Analytic results 
are much harder to obtain and there are finite-size scaling 
problems with numerical methods. Ladder models (dou- 
ble chains) are an interesting intermediate step between 
ID and 2D system. They are easier to treat numerically 
than 2D systems and show new phenomena, which are 
not present in the ID chains. 1-6 Another reason for spe- 
cial interest in ladder systems is the possibility of realiz- 
ing a lattice of weakly coupled ladders in the compounds 
Sr 2 Cu 4 6 and (VO) 2 P 2 7 . 7 ' 8 

A simple but interesting model is the Heisenberg lad- 
& der, consisting of two coupled spin-i Heisenberg chains 
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Here S 8 ' a is the spin operator at site i (i = 1, . . . L) on 
the rung a (a = 1,2) and periodic boundary conditions 
are used along the ladder (^-direction), (see Fig. 1). h 
is an external field in the z-direction. The field h = 0, 
except to calculate numerical derivatives with respect to 
the external field h, and we set gus = 1- The exchange 
constants, J and J', are positive, corresponding to anti- 
ferromagnetic coupling. As the system is translationally 
invariant in the ^-direction the momentum k x is a good 
quantum number. In the ^-direction we use open bound- 
ary condition. This is however equivalent to periodic 



boundary conditions in the ^-direction and a coupling 
J'/2. The momentum along the rungs k y = 0, 7r is there- 
fore also well defined. 

The Heisenberg ladder shows a completely different be- 
havior than the single chain model. While the excitation 
spectrum is gapless (des Cloiseaux - Pearson mode 9 ) for 
the spin S = \ single chain, there exists a spin gap in 
the ladder, 1-3 similar to the Haldane gap in the S = 1 
antiferromagnetic Heisenberg chain. 10 

The Heisenberg ladder and related models, such as the 
t-J ladder 2,4 or the Hubbard ladder 5 and most of the in- 
teresting strongly correlated quantum systems cannot be 
solved analytically. Because of strong interactions mean 
field and perturbation theories often fail to give reliable 
results either. Actually many interesting phenomena in 
these models are of non-perturbative origin. Numerical 
methods giving exact results are thus essential to study 
such systems. 

Four different methods are often used to obtain "ex- 
act" numerical results for strongly correlated systems. 
These are methods without uncontrolled approximations. 
Two of these methods, quantum Monte Carlo (QMC) and 
quantum transfer matrix (QTM) work at finite temper- 
atures. The other two, exact diagonalization (ED) for 
small systems and the density-matrix renormalization- 
group technique (DMRG) are zero-temperature methods. 

The QMC and QTM 11-14 methods are both based on a 
Trotter-Suzuki decomposition of the partition function. 15 
The d-dimensional quantum system is mapped onto a 
(d + l)-dimensional classical one. For quasi ID quantum 
systems, such as chains or ladders, the partition function 
can then be obtained by the QTM. 11 14 This method is 
very powerful. It allows the calculation of the tempera- 
ture dependence of thermodynamic quantities as well as 
correlation lengths for infinite systems. No extrapolation 
is thus necessary for the system size. It does not suffer 
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from the "negative sign problem" of QMC and has much 
higher accuracy. 

We have combined the usual QTM method with 
the look-ahead Lanczos algorithm for non-Hermitian 
matrices 16 to calculate the the extreme eigenvalues and 
eigenvectors of the QTM more efficiently. This new 
method allows us to calculate every thermodynamic 
quantity with higher accuracy. From the numerical point 
of view it is much more efficient then the power methods 
that have usually been used. 

In QMC 13 the partition function is calculated by statis- 
tical sampling of the corresponding classical system, in- 
stead of being evaluated exactly. QMC is very powerful if 
the "negative sign problem" is not severe. It can be used 
in any dimension, on systems with more than hundred 
sites and at lower temperatures than the QTM methods. 
The results are, however, not as accurate as the QTM 
results due to statistical errors from the sampling. These 
errors can be made quite small unless the system inves- 
tigated suffers from the "negative sign problem" . This 
sign problem, which occurs in many frustrated spin sys- 
tems and in 2D fermion systems, often makes simulations 
practically impossible. 

Exact diagonalization by the Lanczos algorithm 17 is a 
very accurate zero-temperature method. It can be used 
to obtain the ground state and the low lying excitation 
spectrum for small systems (of up to about 10 8 states) 
with high accuracy. However the restriction to small sys- 
tems often leads to difficulties with finite size scaling. ED 
can also be used to calculate finite temperature proper- 
ties. But this requires the calculation of a significant 
portion of the energy spectrum or of all energy eigenval- 
ues. The QTM method in contrast needs just a few of 
the extreme eigenvalues of the transfer matrix. 

DMRG 18 is another zero-temperature method. It can 
be used to calculate the ground state and the low ly- 
ing spectrum for larger systems (about 100 sites). This 
method works exceptionally well for one-dimensional 
chains. It can also be applied to higher dimensional sys- 
tems, but there it is harder to obtain accurate results. 

The Heisenberg ladder was studied by Dagotto et al. 2 
and Barnes et al. 3 using exact diagonalization of ladders 
with up to 2 x 12 sites and QMC on systems with up to 
2 x 32 sites. From their finite size results they extrap- 
olated a spin gap of 0.5J for the infinite length ladder 
at the isotropic point J = J'. A calculation by Noack 
et al. 6 using the density matrix renormalization group 
technique 18 (DMRG) gives a spin gap of A « 0.5037J 
and a correlation length of £ = 3.19(1) for J = J'. 

The Heisenberg ladder was also treated in a mean-field 
approximation by Gopalan et al.. 19 They calculate the 
spin gap and the excitation spectrum. The dispersion of 
the spin-triplet excitations agrees well with ED results. 

Using the QTM method we have studied the tempera- 
ture dependence of the correlation length £, the suscep- 
tibility x an d the specific heat C directly for the infinite 
ladder for temperatures down to T ps 0.2J. The spin gap 
and the temperature dependence of £, x, C and of the nu- 



clear spin relaxation rate 1/T\ at low temperatures was 
calculated by combining the QTM with ED results on 
the excitation spectrum. 

II. THE QUANTUM TRANSFER MATRIX 
METHOD 

The QTM method has been widely used to study spin 
models numerically. 11-13,20 The method is based on a 
mapping of the rf-dimensional quantum mechanical sys- 
tem onto a d+l dimensional classical one. For some mod- 
els, e.g. Bethe ansatz solvable models, the partition func- 
tion of the corresponding classical model can also be cal- 
culated analytically by using a transfer matrix method. 
The QTM of a ID spin- 1/2 model is equivalent to the 
diagonal-to-diagonal transfer matrix of the eight-vertex 
model, 21 which can be solved exactly using Bethe ansatz. 
ID models treated analytically include the Heisenberg 
model, 22 the XXZ-model and XYZ-model 23,24 and the 
Hubbard model. 25 

The first step of the QTM is the Trotter-Suzuki de- 
composition of the grand canonical partition function of 
a quantum model. 15 The Hamiltonian is decomposed into 
two parts H = Hi + H2, each of which is easy to diago- 
nalize. A standard choice is the decomposition into two 
sums of commuting terms: 

Hi = H {1 \ Hi = E H(i) > ( 2 ) 

i even i odd 

with [H^ l \ if( J )] = for i,j both even or both odd. 
The two sums Hi and if 2 do not commute in general. 
The simplest decomposition for a chain with only near- 
est neighbor interactions is the so-called "checkerboard 
decomposition". 26 There all terms on odd-numbered 
bonds are collected in Hi and the even-numbered bonds 
into H2 (see Fig. 2a). This is the standard decomposition 
used in most calculations. We have used it in this paper 
for the ID chains. For the ID Heisenberg model the ff W 
are: 

H& = JSi-S i+ i-|(Si + S i+ i). (3) 

A similar decomposition, shown in Fig. 2b, can be used 
for ladder models. For the Heisenberg ladder it is: 

fjM = J ^ Si ta ■ S 8 '_|_i )(I 
a = l,2 
J' 

+ — (Sj',1 • S 8 ' 2 + S 8 '_|_i i • S 8 '_|_i 2) 

~\ E (Si.a + S i + 1 , a ) . (4) 

a = l,2 

Using this decomposition it is possible to approximate 
the partition function in the following way: 
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Z = Ti(e-^ H ) = Ti((UiU 2 ) M ) + 0(Ar 2 ) 
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where /? = l/T denotes the inverse temperature (imagi- 
nary time), M is the Trotter number and Ar = . The 
\ik) are a complete orthonormal system of the states, 
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Note that all the factors in each product commute with 
each other. Since Hi and H2 are chosen to be easy to 
diagonalize the evaluation of the matrix elements (i\Ui \ i'} 
is straightforward. 

The decomposition leads to a systematic error which 
is of order Ar 2 oc M~ 2 . We can extrapolate to Ar — ► 
(M — ► 00) by fitting the results for different Trotter 
numbers M to a polynomial in Ar 2 . 

The above equation (5) can be interpreted as an evolu- 
tion in imaginary time (inverse temperature, also called 
"Trotter" direction) of the state by the "time evo- 
lution" operators Ui and U 2 . Within each time interval 
Ar the operators Ui and and U 2 are each applied once. 
This leads to a graphical representation of the sum on 
a square lattice, where the applications of the operators 
= exp(— Ar_ff( J )) are marked by shaded squares (see 
Fig. 3). The configuration on each time slice corresponds 
to one of the states in the sum (5) for Z. 

The QTM exchanges the space and imaginary time 
direction, problem is reformulated in terms of column- 
to-column transfer matrices Vi and V 2 as shown in Fig. 
3. The partition function can be written similarly to Eq. 
(5) as 

Z = Tr [(V^) 172 ] + 0(Ar 2 ) = Tr(V L/2 ) + 0(Ar 2 ), 

(7) 

where V = V1V2 and L is the length of the chain. Here 
we have used periodic boundary conditions in the space 
direction. Again the transfer matrices are products of 
sparse matrices 



vi= n y(i) > 

1<8<2M; i odd 



(8) 



l<i<2M: i even 



The matrix yW can be calculated quite simply from 
the corresponding matrices Let ci and a 2 denote 

the states on the lower left and right corners of a square 
and let Ti and t 2 denote the states on the upper corners, 
as shown in Fig. 4. Then 
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In order to describe most of the thermodynamic prop- 
erties of the system it is enough to know the extreme 
eigenvalues and eigenvectors of the QTM. This follows 
from an interchangeability theorem, 14,27 which allows us 
to interchange the limit of system size L — ► 00 and the 
limit of Trotter number M —> 00. The free energy density 
(per site or per rung for a single chain or ladder resp.) 
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InZ in the thermodynamic limit is: 



/ = - lim lim — lnTr(V i/2 ) 

= lim In Ai , 

where Ai denotes the largest eigenvalues of V . As we will 
see later the ratio of the two largest eigenvalues deter- 
mines the correlation length of the most dominant fluc- 
tuation: 
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All thermodynamic quantities can be calculated as 
derivatives of the free energy. The magnetic susceptibil- 
ity could, for example, be calculated as a second deriva- 
tive of the free energy density / with respect to the mag- 
netic field h. However numerically it is much better to 
calculate it just as a simple derivative of the magnetiza- 
tion. Indeed it is possible to calculate local quantities, 
such as the magnetization or the internal energy directly 
from the eigenvectors. 

Let us first consider the thermal average of a local 
quantity such as the a-component of the spin S a , the 
particle density n or the energy density. We will call this 
observable we want to calculate A. We can calculate the 
thermal average of this quantity anywhere on the lattice 
due to translational invariance. The effect of the mea- 
surement is to locally change one of the weights: 

(A) L = ±Tr(A 1 V 2 V L ' 2 - 1 ) = ^TriAV^ 2 - 1 ), (11) 

where A = AiV 2 . The matrix Ai is Vi with just the 
matrix V (r > altered: 



n 



(12) 



3<s'<2M; i odd 



where ^l^ 1 ) is the matrix V^ 1 - 1 modified by the measure- 
ment: 



(<t 2 , r 2 |^4 (1 Vi) r i) = ( 0-1,0-2 



Tl,T 2 



(13) 



To simplify this further we can rewrite the trace in 
terms of the right and left eigenvectors \ipf~) resp. {ipf\ 
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of the transfer matrix V . Let us again exchange the limits 
M —> co and L — ► oo. For simplicity we will not write 
the limit lim/vf^oo in the following equations, but it is 
always assumed that this limit is taken. The application 
of the transfer matrix V projects out the eigenvector of 
the largest eigenvalue Ai in the limit L — ► oo: 



lim 

L— ^oc 



Zi(<i>i\v L,4 -Av L ' 4 - 
Ei(<f>i\v L ' 2 \<f>i 



(14) 



Thus local quantities are easy to obtain from the eigen- 
vector corresponding to the largest eigenvalue. The spe- 
cific heat C can now be calculated by a numerical deriva- 
tive of the internal energy: 



C 



d(H) 
dT ' 



(15) 



The magnetic susceptibility x can be calculated as a nu- 
merical derivative of the magnetization with respect to 
the external field h: 



d(S z 



dh 



(16) 



h = 



Similarly we can calculate correlation functions, such 
as spin correlations. Let us calculate the correlations of 
the fluctuations of such a quantity around its mean value, 
Ai = Ai — (A), between sites i and i + d. In the limit 
L — ► oo this is: 



{AiAi +d ) 
for odd d, and 
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for even d. These correlations are simple to calculate for 
short and intermediate ranges d. Often more interesting, 
and much simpler to calculate, is the correlation length, 
defined as 
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As we want to take the limit d —> oo it is sufficient if we 
consider the case of even d. In the limit d —> oo formula 
(18) becomes 



lim (AiA i+d ) 

d^-oo 



lim 

d— ^oo 



lim 

d— ^oo 



(^\A\^)(^\A\^) 
(^f)AiA a 

(tf\A\1>2){rt\A\1>?) 
(^|^f)AiA a 



Ai 



d/2 



■ exp 



ikd 



(20) 



A„ is the largest eigenvalue with nonzero overlap 
(i>i\A\i>a)(i>a\A\i>?-}. If the state A\ipf) is in the same 
invariant subspace as (e.g. if A = S z ), then it is 

usually the second largest eigenvalue in this subspace, 
otherwise (e.g. if A = S x or A = S y ) it is usually the 
largest eigenvalue in the invariant subspace that contains 
A4\ipf~) . In Eq. (20) it was assumed that there is only one 
eigenvalue with absolute value |A„|. The generalization 
of the above formula to the case of multiple eigenvalues 
with the same absolute value (e.g. a complex conjugate 
pair) is straightforward. 
The correlation length £ is 
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and the wave vector of the most dominant fluctuation k 
can be calculated from the phase of A„: 



7 r 1 (A a 
k = lim — arg — 

m^oo 2 VAi 



+ ri7r (n = or 1). (22) 



The ambiguity arises because the transfer matrix in this 
formulation propagates over two sites and cannot dis- 
tinguish between k and k + ir. It can be resolved by 
comparing the correlations for odd and even d. 



III. THE LOOK-AHEAD LANCZOS ALGORITHM 

The numerical problem in the QTM method is the cal- 
culation of the extreme eigenvalues and the correspond- 
ing eigenvectors of the transfer matrix V . This is a sim- 
ilar problem as in exact diagonalization (ED). In ED we 
want to calculate the lowest eigenvalues and eigenvectors 
of the Hamiltonian. ED is restricted to small system 
sizes, as we have to store three vectors of the Hilbert 
space in the main memory of the computer. In the QTM 
method we have exchanged the space direction with the 
imaginary time direction. The length of the chain can 
now be made as large as one wishes. The price we have 
to pay is that we have to store the vectors of possible 
states in the imaginary time direction. We are restricted 
to a small number of time slices and thus to the high and 
intermediate temperature regime. 

The main problem is that, while both V\ and V 2 are 
hermitian, their product is no longer hermitian, since the 
two matrices do not commute. Until recently there was 
no efficient way to calculate eigenvalues and eigenvectors 
of non-hermitian matrices, since the usual Lanczos algo- 
rithm is numerically unstable for non-hermitian matrices, 
and usually does not converge. Therefore the eigenvalues 
and eigenvectors were calculated using power methods. 
Recently however a variant of the Lanczos algorithm, the 
look-ahead Lanczos algorithm was developed. 16 This is 
almost always numerically stable and convergent. Very 
rare exceptions, so-called "incurable breakdowns", can 
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usually be circumvented by using different starting vec- 
tors. We have never encountered such an incurable break- 
down in our calculations. 

The Lanczos algorithm 16,17 is an iterative method to 
tridiagonalize a matrix V . The extreme eigenvalues of the 
recursively generated tridiagonal matrix converge very 
rapidly to the eigenvalues of the original matrix. As the 
matrix V is needed only in form of matrix-vector prod- 
ucts Vv the Lanczos algorithm is ideally suited to cal- 
culate the extreme eigenvalues and eigenvectors of large, 
sparse matrices. 

The Lanczos algorithm recursively generates the tridi- 
agonal matrix and two sets of vectors {i> 8 } and {u> 8 } 
(i = . . .N — 1) starting from the vectors vo and wq. 
These basis vectors span the N-ih Krylov subspace of V 
and 0: 

span({i; 8 }) = span{i; , Vv , . . . , V N ~ 1 v }, (23a) 
span({w 8 }) = span{w , V^w , • • • , (V^) N ~ 1 w }, (23b) 

and they are biorthogonal: 

(v i ,w j ) = b ij . (24) 

The Lanczos algorithm terminates regularly when an 
invariant subspace of V or has been found and vjy = 
or u>jv = 0. For non-Hermitian matrices a breakdown 
occurs, when the the vectors vn and wjy are orthogonal. 
Then vn ^ and ^ 0, but (vn,wn) = and the 
normalization (Eq. (24)) cannot be fulfilled. In finite 
precision arithmetic there can also be near-breakdowns, 
when i>jv and wjy are nearly orthogonal and the algorithm 
becomes numerically unstable. 

The Lanczos algorithm is most often used for Hermi- 
tian matrices, where such breakdowns cannot occur. If 
we choose vo = u>o then we have = for all i since 
V = Vt. The normalization Eq. (24) then becomes sim- 
ply: 

(v N ,w N ) = (v N ,v N ) = \\v N \\ 2 . (25) 

This is zero only in the case of regular termination, where 
v _/y = wjv = 0. 

The look-ahead Lanczos algorithm relaxes the condi- 
tion of tridiagonalizing the matrix. As long as there are 
no breakdowns or near-breakdowns it is equivalent to the 
usual Lanczos algorithm. If a breakdown would occur 
in calculating vn and it tries to skip over that it- 
eration. The simple three-term recurrence relations of 
the standard Lanczos algorithm are then replaced by 
more complex relations including not only the vectors 
i>iv_2, i>jv-i, Vvn-i and wjv-2, wn-i, V^w^-i but also 
V 2 v N _ 1 ,. . . , V'vN-uiVlfwN-u (V^'wn.u . . .. I 
is the look-ahead length. The look-ahead Lanczos al- 
gorithm then generates a block-tridiagonal matrix with 
blocks of size / instead of a tridiagonal one. Usually a 
look-ahead of / = 2 or 3 is sufficient except in rare cases. 
In extremely rare cases we would encounter breakdowns 
with any number of look-ahead steps. This case is called 



an incurable breakdown. For details we refer to the origi- 
nal literature. 16 An implementation of the eigenvalue al- 
gorithm is available in electronic form. 28 

The look-ahead Lanczos algorithm allows us to calcu- 
late the extreme eigenvalues of the QTM very efficiently 
and with high accuracy. We need much less iterations 
compared to the power method. We found that the look- 
ahead Lanczos algorithm often converges in just a few 
dozen iterations. 

Another advantage is that the eigenvectors can be cal- 
culated without any problems by the Lanczos algorithm. 
This allows us to calculate quantities such as the inter- 
nal energy or the magnetization directly via Eq. (14). 
These results are more accurate than the calculation as 
numerical derivatives of the free energy. 

We have compared the algorithm to exact results for 
the 1-D XY and Heisenberg models. 24 We found that our 
results are very accurate down to quite low temperatures 
(T pa 0.1 J) for results extrapolated from M = 1 . . . 10. 

IV. RESULTS FOR THE HEISENBERG LADDER 
A. Quantum transfer matrix results 

As an application of the new algorithm we have stud- 
ied the Heisenberg ladder. Specifically we have calculated 
the correlation length £, the specific heat C and the mag- 
netic susceptibility x as a function of the temperature T. 

In Fig. 5 we show the susceptibility per spin x as a 
function of the temperature for different values of J/J'. 
At high temperatures the results agree well with a third- 
order high temperature expansion: 

X (T) = {T- 1 - |(J + \J')T- 2 + £ JJ'T~ 3 . (26) 

At low temperatures we observe an exponential drop of 
the susceptibility, caused by the gap in the spin excitation 
spectrum. This drop is steeper for smaller values of J/ J', 
indicating that the gap A/J' decreases with increasing J. 
The spin gap will be studied in more details in Sec. IV B. 

The spin gap also leads to an exponential drop of the 
specific heat, as shown in Fig. 6. At high temperatures 
there is again good agreement with the high-temperature 
expansion. The free energy per site is 

/(T) = -Tln2-^(j 2 + ijJ')T- 1 , (27) 

and the specific heat 

C(T) = ^(J 2 + \JJ')T- 2 . (28) 

In Fig. 7 we show the temperature dependence of the 
correlation length £ for the Heisenberg chain and the 
Heisenberg ladder, calculated by the QTM. The wave 
vector of the dominant correlation is k = [tt, 7r) for 
the ladder, which corresponds to antiferromagnetic cor- 
relations. In the high temperature limit the correlation 
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length is similar in both models. With decreasing tem- 
perature the correlation length becomes longer for the 
ladder. This is because antiferromagnetic correlations are 
enhanced faster in the ladder due to the larger number 
of nearest neighbor sites. At low temperatures the cor- 
relation length saturates to a finite value, £ ps 3-4, which 
agrees with £ ps 3.19, determined by the DMRG calcula- 
tion for zero temperature. 6 This finite correlation length 
corresponds to a finite spin gap via the relation, £ ~ a/A, 
where a is a constant of the order of a characteristic spin 
velocity. In the gapless single chain, on the other hand, 
the correlation length diverges like £ ~ v s /(ttT) (v s : the 
spin velocity) as predicted by conformal field theory. 29,30 



With increasing J the collective one-"magnon" branch 
crosses into the two-"magnon" continuum (see Fig. 8c). 
The exact diagonalization 3 and mean-field 19 results indi- 
cate that even then the spectrum can still be described 
by the above picture. 

Using these results on the excitation spectrum we can 
calculate the low temperature thermodynamics of the 
Heisenberg ladder. First we start from the simple limit 
J = 0. Each rung can be either in the singlet state or in 
one of the three triplet states. We obtain for the partition 
function of the ladder of length L: 



Z = (l + e -/3(A + ft) + fi -/3A + fi -/3(A-ft))£ 

= [l + (l + 2cosh(/3/ l ))e-' 3A ] 1 ', 



B. Spin gap and low-temperature thermodynamics and for thg magnetic suscep tibility 



To calculate the spin gap and the thermodynamic 
quantities at low temperatures we start from the limit 
J/J' — ► 0, where a simple description of the whole ex- 
citation spectrum is available. In that limit, each eigen- 
function of the total system can be written as a direct 
product of one-rung states, which are either spin singlets 
or one of the triplets (a = —1, 0, 1), and the ground state 
is that with all singlets. Accordingly, each eigenenergy 
is given by J'N, where N is the number of triplet rungs, 
measured from the ground state energy — | J'L, and the 
energy spectrum shows a tower structure consisting of 
equidistant multiplets with separation J'. Each multiplet 
is labeled by the number of triplet rungs, N . The first 
excited multiplet consists of the states with one triplet 
rung and therefore belongs to the sector of Stot = 1; 
and its multiplicity is ?>L. In general, the 7V th -multiplet, 
which consists of the states with N triplet rungs, has the 

multiplicity, g(L,N) = 3^ (jv)> where the first factor 

3^ comes from the spin part. 

A small but finite value of J lifts the degeneracy of 
these states. A schematic picture of the energy levels is 
shown in Fig. 8. The one-triplet excitations then form a 
three-fold degenerate band of collective excitations with 
dispersion e k = J' + J cos k and z-component of spin a = 
— 1, 0, 1, where we set the ground state energy Eq.s. = 0. 
The minimum of this band is at a momentum k x = 7r 
along the ladder. The momentum along the rung is k y = 
7r. We will call these excitations "magnons" although 
there is no magnetic long range order in the ground state, 
order . To second order in perturbation theory the gap 
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j + \jyj>. 



(29) 



The higher excited states form a continuum of excited 
states, with k y = and a minimum at k = 0. They can 
be viewed as two-"magnon" states. In low-order pertur- 
bation theory the magnon-magnon interaction is repul- 
sive. The minimum of the continuum is thus at energies 
slightly larger than twice the gap 2 A. 
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which drops like /3e~ l3A for low temperatures. 

If J is nonzero we have to take into account the disper- 
sion of the spin excitations. We assume the magnon exci- 
tations to have a dispersion +crh, where a = — 1, 0, 1 is 
the z-component of the spin. In the limit T —> the inter- 
actions between these magnons become negligible since 
the magnon density goes to zero due to the gap. 

The "magnons" are boson-like in the sense that one 
can excite more than one excitation with the same quan- 
tum numbers, the wave number k and the spin a, but 
they are not real bosons, since the Hilbert space is re- 
stricted. One cannot excite two "magnons" at the same 
rung, which might be described by a hard-core repulsion 
in the boson representation. This was first pointed out 
by Dyson for real magnons in a ferromagnetic state, usu- 
ally referred to as kinematical interactions. 32 The kine- 
matical interactions become important with increasing 
temperature, and are essential to get a correct tempera- 
ture dependence. Otherwise, for example, in the limit of 
T —> oo the number of bosons at each k would diverge 
and we would not obtain the correct entropy for T —> oo. 

At low enough temperatures, T <C A, the magnon 
density is very low and it is sufficient to include up to 
one magnon for each k and a. Therefore, both residual 
magnon-magnon interactions and the kinematical inter- 
actions are negligible. The free energy per site in that 
limit is 



/l = -^ [1 + 2c ° sh ^ ] E' 



2L(] 
1 



2/3 



[l + 2cosh(/3/i)]z(/3), 



(32) 



where we have replaced the sum by the integral 

I /»7T POO 

z(/3) = — dke- 13 '" = Acp{c)e- I}e , (33a) 
2t J.jt Jo 
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which is a Laplace transform of the magnon density of 
states, p{e). The susceptibility then is 



xi CO 



5 2 /i 



dh 2 



(34) 



h = 



For a simple form = A + a | \k \ — ir\ n we can perform 
the integration 



z(p) 



r(i) 



(a/3) 



-1/n -/3A 



(35) 



where we extended the integration over A; to infinity. At 
low temperatures the magnetic susceptibility then is 



x ( T ) = £iAl a -l/n T -l + l/n g — A/T ^ 



(36) 



If we replace the magnon band by the quadratic approx- 
imation (n = 2) we get 



X(2)(T) 



1 



-A/T 



Similarly we can calculate the specific heat as 



(37) 



C(n)(T)= — - 



^\ l/n /rp\2—l/n 



A 



-A/T 



r(i) + 2r(i + i)| + r(2 + i) (J 



(38) 



In the quadratic approximation 

C\2)(T) 



3 / A \ / \ 



4 V 7ra 



A 



T 3 /T 



A 4 V A 



-A/T 



(39) 



The low-temperature result Eq. (36) motivates a first 
estimate of the gap based on the logarithmic derivative 
— dg 1 ^ ■ This derivative is A — (1 — \/n)T at low tem- 
peratures for a susceptibility (36). This derivative goes 
to zero on the other hand if the susceptibility is nonzero 
for T = or vanishes following a power law. 

In Fig. 9 we show — for some gapless systems, 
the ID Heisenberg and XY-models, to compare with the 
Heisenberg ladder. This plot clearly shows the existence 
of a spin gap for the ladder. 

The size of the gap however is not easy to estimate from 
the data. For J = J' we can reach only temperatures 
T/J' pa 0.2, which is below the gap but not yet in the 
asymptotic region. 

To determine the size of the spin gap more precisely, 
we need a fitting function which describes the whole tem- 
perature range. This function should give correct results 
in both low and high temperature limits. To get a cor- 
rect high temperature limit, one has to take into account 



the kinematic interactions, as discussed before. We have 
found a simple way of including kinematical interactions 
in the thermodynamics based on reasonable physical ar- 
guments. Our formula not only gives correct low and 
high temperature limits, but the overall agreement also 
turns out to be nice. 

In our new formula, the grand partition function is 
calculated as follows. The main problem of the boson 
description is that as the number of triplet rungs, N, 
increases, the number of the corresponding boson basis 

states diverges like gg(L,N) = ( 3L+ n~ 1 ^) > wnl l e the 

correct dimension is g(L,N) = 3^ (jy)' Therefore, the 
basic idea is to reweight the 7V-magnon part in the parti- 
tion sum, [g(L, N)/ ' gs(L, N ]\Z\- Joson (N -ra&gnon) , so that 
each multiplet contributes the correct entropy. 33 This can 
be done with slight modification of the boson partition 
function, 



Z' C (L,N)= J2 



exp 



{kj,aj} 



N 



(40) 



This corresponds to a sum neglecting the undistinguisha- 
bility of bosons. Aside from the global factor, the differ- 
ence from the original boson sum are the terms in which 
two or more bosons have the same quantum numbers, 
(k, a). However, the number of these terms is smaller by 
at least order one w. r. t. L and we neglect those correc- 
tions. Since there are (3L) N terms in Z' c , the reweighting 
should work as follows: 



2 = E 



N=0 
L 



9 -N±Z'c(L,N) 

(3L) N CK ' ; 



N=0 



J2[ N ) L ~ N £ (l + Scosh^re-^,^-, 



ki,...,kp. 



l + (l + 2cosh(/3/ l ))|^ 



(41) 



There are 4 L terms in total in the above partition sum, 
giving the correct total entropy, since we reweighted to 
get the correct number of excitations. Note that here 
(i) we assume that all excitations could be described as 
multi-magnon excitations and (ii) all residual magnon- 
magnon interactions are neglected. The assumption (i) 
is obviously correct in the limit of J/ J' — ► and there is 
no indication of a breakdown of the arguments of analytic 
continuation w. r. t. J: e.g., the spin gap is always finite 
as far as J is nonzero. 

The free energy per site is 



/ 



1 



In [l + (l + 2cosh(/3/i))z(/3)] 



(42) 



where we have taken the limit L — ► oo and again replaced 
the sum over k by an integral. 
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This partition functions gives a susceptibility 



X = P 



z(p) 



1 + 3z(/?) : 



(43) 



which is correct in both limits T —> and T —> oo. For 
very low temperatures we recover the result of the low- 
temperature approximation (34). For high temperatures 
we obtain the correct Curie law \k, 

We will now try to fit the QTM results to the above 
model. The function z(JS) depends on the dispersion ej, 
we use. First we discuss the small J/J' region. As the 
correction term in the Trotter-Suzuki decomposition is 
of order (3 3 J 3 / M 2 we can reach quite low temperatures 
T/J' when J <C J' ■ For J = 0.1J' we can reach temper- 
atures below T/J' pa 0.04. These temperatures are low 
enough to see the asymptotic behavior. In that limit the 
dispersion is of the form 



J' + J cos k = A + 2J cos 2 (*/2), 



(44) 



with A = J — J' . For J/ J' = 0.1 we can reach quite low 
temperatures and a fit to the above Eq. (43) using the 
dispersion (44) is excellent. The resulting fit of x(T) is 
shown in Fig. 10. A least square fit gives a gap of A = 
0.909, which is in excellent agreement with the second 
order perturbation result A = 0.905 (Eq. (29)). 

At J = J' the gap is harder to estimate. This is caused 
by two facts. First we cannot reach as low temperatures 
as in the small J region, as the Trotter "time" step /3J/M 
is now larger. The lowest temperatures we can reach are 
about T/J pa 0.2. Additionally we do not know the exact 
shape of the dispersion. 

We can guess the form of the dispersion from exact 
diagonalization data 3 and mean-field calculations. 19 The 
dispersions obtained by both the perturbation result Eq. 
(44) and the mean-field are quadratic close to the min- 
imum at k = tt. At larger \k — tt\ exact diagonalization 
and mean-field results indicate a more linear behavior. 
We have used several functional forms for the dispersion. 
Good fits were obtained by the following dispersions: 



e (1) 

t k 
J 2 ) 



^A 2 + 4Aa(l + cosfc) 2 
^A 2 + 2Aa(k - tt) 2 , 
A + a(\k\ - tt) 2 



A- 



e (3) . 
t k ' 

4 4) = A + c||fc| 



if ||*| 



7T < 



+ c||*| 
-tt|. 



7r| otherwise, 



(45a) 
(45b) 

(45c) 
(45d) 



The dispersion (45a) is the functional form obtained by 
the mean-field calculation. 19 

In Table I we show the gap, the curvature a = 

f afc 2 an d the other fitting parameters obtained by 

a least square fit of the QTM results for In x to the above 
dispersions (45). 

The discrepancies between the fits arise because, al- 
though we can simulate at temperatures below the gap 
T ps 0.4 A, we are not really in the low-temperature 



regime where the interactions between the excitations 
and the exact shape of the dispersion become unimpor- 
tant. This can be seen best in the plot of — in Fig. 
9. In Fig. 11 we show the fit of the susceptibility for the 
dispersion (45a). The susceptibilities obtained using the 
other dispersions differ only slightly. 

The dispersion e^ 4 " 1 is not realistic, as it is not 
quadratic, but linear close to the minimum at k = tt. 
It underestimates the gap, since the density of states is 
too small near the minimum. For the same reason we 

(3) 

believe that the dispersion e k underestimates the cor- 
rect gap. Similarly a dispersion that is too flat close to 
the minimum overestimates the density of states there 
and thus also the gap. We estimate the gap to be in the 
range 0.45J < A < 0.5J, which is in agreement with the 
exact diagonalization 3 and DMRG results. 6 



C. Nuclear spin relaxation rate 

Another quantity of interest is the nuclear spin relax- 
ation rate \/T\. It can be written in terms of the dy- 
namical susceptibility perpendicular to the field: 34 



±- = 2 7 2 Tj2\A 



iX"(q,wo) 



LO 



(46) 



where \A C 



(2L) 1 is the form fact' 



or, 7 is 



the 



clear gyro-magnetic ratio, and uio the nuclear resonance 
frequency, which is a very small energy scale, typically of 
the order of mK. The main problem here is the calcula- 
tion of the imaginary part of the susceptibility x"(q, wo). 
This can be related to the dynamical structure factor by 
the fluctuation-dissipation theorem: 



X'_[(q,w ) = <S|_(q,Wo)(l 



)«Sj.(q,wo)/?wo, (47) 



where we have used the fact that u>o ps 3mK <C T. As 
the Hamiltonian of the system is invariant under spin 
rotations and there is no long range order present the 
susceptibilities in all directions are equal: xj_ = Xzz and 



<Sj_(q, ojq) = S z (q,u) ) 



E 



\(m\S z Jn)\ 2 6(E m - E n 



(48) 

■w )e-(> E "/Z, 



where |m), \n) are complete sets of eigenstates with en- 
ergy E m and E n respectively, and 



Si 



1 



r qz 



(49) 



Which states contribute to \/T\ at low temperatures? 
The dominant fluctuations in the ground state are anti- 
ferromagnetic, leading to a maximum in the equal-time 
spin structure factor at q = (tt,tt). However these dom- 
inant antiferromagnetic fluctuations (m|5'q|G.5'.) of the 
ground state near q = (tt, tt) do not contribute since they 
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have an energy gap of E m — Eq.s. > A ^ wo- The only 
relevant contributions arise from fluctuations of the ex- 
cited states with small momentum transfer q. 

At low temperatures we assume the one-magnon states 
to be independent. We restrict the sum over n to the 
independent one-magnon states \k,a) with momentum 
k and z-component of spin a = — 1,0,1. As wo < A 
and momentum is conserved only the states \k + q x ,a), 
contribute to the sum over m: 



S z (q,ui ) x, ^2\(k + q x ,a\S^\k,a)\ 2 6 qyt0 



x8(e k+qx - e k - oj )e 



(50) 



where \k, a) is a one-magnon state with momentum k and 
z-component of spin a = —1,0, 1. From the excitation 
spectrum it is obvious that only for q y = and q x ps 
or q x pa —2k we have a nonvanishing S z (q,uio). Using a 
second order Taylor expansion of the dispersion we can 
write S z (q, ujq) in terms of (^-functions of q: 



S z (q, ui ) & ^2\(k + q x , a\Sq\k, a)\ 2 6 qyt0 e- 



I3e k 



k,a 



6(q x ) + 8(q x + 2k) 



V (k)\ 1-2^^, 



(51) 



dk 



where we have set loq —> in the (^-functions. v(k) — ^ L 
is the group velocity. The matrix elements are: 

1 , 



\(k,a\S( 



0,0) 



k, a) 



2L 



(52) 



We estimated the matrix element |(— k, u\S z _ 2k ol^> <7 )| 2 
by exact diagonalization on finite ladders of up to 10 
rungs. For J = J' is is nearly constant for tt/2 < \k\ < tt: 



\{-k,*\S(_ 2kA) \k,*)\ 2 x0.5 



1 

2L C 



(53) 



At low temperatures the main contributions arise from 
k pa 7r, where q ?a (0, 0). We replace the matrix elements 
\A q \ 2 by its value \A q \ 2 x A 2 = A 2 /2L at q = (0,0). 
Replacing sums by integrals we get in the quadratic ap- 
proximation for the dispersion in the temperature range 
wo < T < A: 



J_ _ 2j 2 A 2 



4tt ^ 

q y = 0,w 

2 A'l 



dk- 



dq x S±(q,ui ) 
3e" £ */ T 



8tt 2 J UA, 2a^(7r - k) 2 + LO /a 



16a7r 2 



\2T) 



(54) 

(55) 
(56) 



where A'o is the modified Bessel function of second kind. 
In the temperature regime where our approximation is 
valid we can expand A'o(fy) ~ — C + In 4 — ln(uio/T) ps 



0.80908 - ln(w /T). C « 0.577216 is Euler's constant. 
Thus finally we have for the nuclear spin relaxation rate 



Zy 2 A 2 
16a7r 2 



-A/T 



(0.80908 - ln^o/T)), 



(57) 



in the temperature range u>o 3mK <C T <C A. The 
main feature is the exponential drop with temperature 
caused by the gap. In addition there is a logarithmic 
divergence in u>o caused by the van Hove singularity at 
the band minimum in the density of states of spin exci- 
tations. Although the equal time spin correlations have 
a maximum at q = [tt, 7r) these fluctuations do not con- 
tribute since they have a large energy gap. As the main 
contribution to the nuclear spin relaxation rate comes 
from q ps (0, 0), and not from q ?a [tt, 7r) we expect the 
temperature dependence to be similar for Cu and O sites 
in a copper-oxide ladder. This differs from the case of 
copper-oxide planes, where there is a marked difference 
in the temperature dependence, because there are low 
energy fluctuations around q = [tt, 7r) that contribute to 
l/Ti at Cu sites but not at O sites. 35 



V. CONCLUSIONS 

We have developed an improved version of the quan- 
tum transfer matrix algorithm. Quantum transfer matrix 
methods (QTM) do not suffer from the sign problem of 
quantum Monte Carlo. Therefore they are ideal to inves- 
tigate models where the sign problem is severe. Exam- 
ples include frustrated spin systems or fermionic ladder 
models, like the t-J ladder. 

We have combined the QTM method with the look- 
ahead Lanczos algorithm to calculate the extreme eigen- 
values and eigenvectors of the QTM with very high ac- 
curacy. The algorithm converges much faster than usual 
power methods. The calculation of the eigenvectors of 
the transfer matrix with high precision by the look-ahead 
Lanczos algorithm allows a direct calculation of the mag- 
netization, internal energy, magnetic susceptibility, spe- 
cific heat and similar quantities for an infinite length sys- 
tem. 

In this paper we have reported on the thermodynam- 
ics of the Heisenberg ladder. The QTM method by it- 
self is restricted to high and intermediate temperatures 
(T > 0.2J). By combining the QTM method with exact 
diagonalization results for the low-lying excitation spec- 
trum we are able to calculate the temperature depen- 
dence of the specific heat, magnetic susceptibility and 
the correlation length for the entire temperature range. 
This also allows an estimation of the spin gap. Finally 
we have calculated the temperature and frequency de- 
pendence of the nuclear spin relaxation rate l/Ti. 

In the final stages of the preparation of the manuscript 
we learned about a preprint by Barnes and Riera, in 
which they report on a calculation of the tempera- 
ture dependence of the magnetic susceptibility by exact 
diagonalization. 36 



9 



An interesting question arising here is, what happens 
to the spin gap upon doping of holes (t- J-ladder model). 
This is currently being investigated 2,4,5 
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FIG. 4. Rotation of the transfer matrix: The matrix 
propagates a state along the imaginary time direction. 
propagates along the space direction. 



FIG. 1. Diagram of the Heisenberg ladder with two legs in 
the x direction and L rungs in the y direction. The coupling 
along the legs is J and the along the rungs J'. 



FIG. 2. Examples of decompositions used in the Trot- 
ter-Suzuki decomposition, (a) the checkerboard decompo- 
sition, the simplest decomposition for ID chains. (b) a 
"checkerboard" decomposition for ladder models. 
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FIG. 3. Graphical representation of the Trotter-Suzuki de- 
composition of a one dimensional quantum chain using the 
checkerboard decomposition. Also shown is the formulation 
in terms of the usual row to row transfer matrices U and in 
terms of column to column transfer matrices V. The matrices 
that are altered for measurements are indicated by a lighter 
shading and are labeled. Refer to the text for details. 
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FIG. 5. Temperature dependence of the magnetic sus- 
ceptibility of the Heisenberg ladder for different values of 
J I J' = 1,0.5,0.2 and 0.1. (a) a logarithmic plot of x as a 
function of the inverse temperature /?. (b) \ as a function of 
the temperature T. 
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FIG. 6. Temperature dependence of the specific heat of the 
Heisenberg ladder for J' = J . ^ 
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FIG. 7. Correlation length of the Heisenberg chain and lad- 
der as a function of temperature for J' = J. In the gapless 
Heisenberg chain £ diverges for T — > oo, while it remains finite 
for the Heisenberg ladder which exhibits a spin gap. 
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FIG. 8. (a) Evolution of the energy levels when the interac- 
tion J is turned on. (b) Qualitative picture of the dispersion 
at small J/ J'; (c) at J = J'. 



12 



1.0 



0.8 - 



0.6 - 



J/J'=0.1 



t QTM results 
fit 1 

fit 2 

fit 3 

fit 4 



0.4 



0-2 - Heisenberg chain 



J/J'=1 



0.0 



-0.2 



0.0 



XY-chain 

0.2 0.4 



0.6 



1.0 



FIG. 9. Temperature dependence of the logarithmic deriva- 
tive of the magnetic susceptibility with respect to the in- 
verse temperature fi. Shown is the magnetic susceptibility 
for the gapless ID XY and Heisenberg chain and for the gap- 
ful Heisenberg ladder with J/ J' = 1 and J/ J' = 0.1. Also 
included is the fit which is described in the text. For J/ J' = 1 
the four different fits are shown. The temperature is in units 
of J for the single chains, and in units of J' for the ladder. 
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FIG. 11. Temperature dependence of the magnetic suscep- 
tibility of the Heisenberg ladder with J' = J. The solid line 
is the fit discussed in the text, (a) x as a function of the 
temperature T; (b) a logarithmic plot of x as a function of 
the inverse temperature fi. 
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FIG. 10. Temperature dependence of the magnetic suscep- 
tibility of the Heisenberg ladder with J/ J' = 0.1. (a) x as a 
function of the temperature T; (b) a logarithmic plot of \ as 
a function of the inverse temperature fi . 
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